BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) u₀ = prnbng.u₀ p = prnbng.p tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, u₀, tspan, p) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, u₀, tspan, p, jac=true)
Parsing parameters...done
Adding parameters...done
Parsing species...done
Adding species...done
Creating ModelingToolkit versions of species and parameters...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, u₀, tspan, p, jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 392s / 99.0% 71.2GiB / 99.5%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
ODEProb DenseJac 1 264s 68.1% 264s 46.8GiB 66.1% 46.8G
iB
ODEProb SparseJac 1 73.0s 18.8% 73.0s 16.7GiB 23.6% 16.7G
iB
ODEProb No Jac 1 34.6s 8.90% 34.6s 5.07GiB 7.16% 5.07G
iB
Parse Network 1 10.4s 2.67% 10.4s 1.22GiB 1.73% 1.22G
iB
Create ODESys 1 5.84s 1.51% 5.84s 1.05GiB 1.48% 1.05G
iB
──────────────────────────────────────────────────────────────────────────
──
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show numparams(rn) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 numparams(rn) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = copy(u₀) du = similar(u) @timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.) @timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.) densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.)
We also time the ODE rhs function with BenchmarkTools as it is more accurate given how fast evaluating f is:
@btime oprob.f($du,$u,$p,0.)
23.810 μs (3 allocations: 512 bytes)
Now we time the Jacobian functions, including compilation time in the first evaluations
J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 824s / 97.5% 85.0GiB / 99.2%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
SparseJac Eval1 1 339s 42.2% 339s 7.20GiB 8.53% 7.20G
iB
ODEProb DenseJac 1 264s 32.9% 264s 46.8GiB 55.5% 46.8G
iB
ODE rhs Eval1 1 74.0s 9.22% 74.0s 4.35GiB 5.15% 4.35G
iB
ODEProb SparseJac 1 73.0s 9.08% 73.0s 16.7GiB 19.8% 16.7G
iB
ODEProb No Jac 1 34.6s 4.30% 34.6s 5.07GiB 6.01% 5.07G
iB
Parse Network 1 10.4s 1.29% 10.4s 1.22GiB 1.45% 1.22G
iB
Create ODESys 1 5.84s 0.73% 5.84s 1.05GiB 1.25% 1.05G
iB
DenseJac Eval1 1 2.41s 0.30% 2.41s 1.96GiB 2.32% 1.96G
iB
SparseJac Eval2 1 62.8μs 0.00% 62.8μs 912B 0.00% 91
2B
ODE rhs Eval2 1 55.2μs 0.00% 55.2μs 688B 0.00% 68
8B
──────────────────────────────────────────────────────────────────────────
──
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)
